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Abstract 

Nose’s pioneering 1984 work inspired a variety of time-reversible deterministic thermostats. 
Though several groups have developed sucessful doubly-thermostated models, single-thermostat 
models have failed to generate Gibbs’ canonical distribution for the one-dimensional harmonic 
oscillator. A 2001 doubly-thermostated model, claimed to be ergodic, has a singly-thermostated 
version. Though neither of these models is ergodic this work has suggested a successful route 
toward singly-thermostated ergodicity. We illustrate both ergodicity and its lack for these models 
using phase-space cross sections and Lyapunov instability as diagnostic tools. 
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I. SINGLE-VARIABLE THERMOSTATS AND GAUSSIAN ERGODICITY 


In 1984 Hoover explored the application of the NosAHoover version^ of Nose’s canonical 
motion equations^>^ to a harmonic oscillator at thermal equilibrium with coordinate q , 
momentum p , temperature T , and thermostat variable ( : 

{q = p ; p= -q-Cp ; ( = [ p‘^ - T ]/r^ } [ NH ] . 

Posch, Hoover, and Vesely found that this model partitions the {q,p,C) phase space into 
many separate toroidal regions embedded in a chaotic sea4. The complexity and the stiffness 
of the solutions increase rapidly as the thermostat response time r is reduced. In addition to 
equilibrium applications analogous motion equations can be used to thermostat irreversible 
nonequilibrium simulations such as steady shear and heat flows. The harmonic oscillator 
can generate steady-state heat flow problems if the temperature varies in space^*^ : 

1 — e<T = T{q) = 1 + etanh(g) < 1 + e . 

Here e is the maximum value of the temperature gradient, {dT/dq), to which the oscillator 
is exposed. It can be viewed as the strength of nonlinearity, and depending on its value, one 
can move from the equilibrium regime (where e = 0) to the nonequilibrium regime (where 
e > 0). 

Somewhat paradoxically, the NosAHoover motion equations as well as all the others we 
consider here are time-reversible, even away from equilibrium. That is, any time-ordered 
sequence of {q,p,C) points can be reversed either by [ 1 ] changing the sign of dt in the 
integrator, or [ 2 ] by changing the signs of the {p, () variables. The harmonic oscillator 
equations also have mirror symmetry. Changing the signs (+g, +p) <—)■ (—q, —p) gives an 
additional pairing of solutions. 

Apart from being time-reversible, a good thermostat must result in ergodic dynamics. 
Ergodicity of the dynamics connects dynamical averages with corresponding Boltzmann- 
Gibbs phase averages. In describing the results of the present work, we have used Ehrenfests’ 
idea of “quasiergodicity”, where the dynamics eventually comes arbitrarily close to each 
feasible point, interchangeably with “ergodicity”. 

For the equilibrium NosAHoover harmonic oscillator, the Gaussian distribution is the 
stationary solution of the Liouville’s phase-space continuity equation: 

V = r = {q,p,C) —^ {.df/dt) = • {fv) =0 —^ 
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On the other hand, nnmerical work gives two kinds of solntions, either quasi-periodic tori, 
with all Lyapnnov exponents being zero or a single chaotic, Lyapnnov-nnstable sea. The 
global dynamics, therefore, either remains conhned within the tori, or occupies the chaotic 
sea separated from the tori, depending upon the initial conditions, the temperature T , 
and the response time r. In other words the presence of two sets of global maximal Lya¬ 
punov exponent - one positive and another zero, indicates that a trajectory starting from 
an arbitrary initial condition is unable to explore the neighbourhood of the entire feasible 
phase-space. As a result, the phase-space gets partitioned into at least two noncommunicat¬ 
ing regions, violating the metric indecomposibility of the phase space - the necessary and 
sufficient condition for ergodic dynamics according to Birkhoff’s theorem. Thus the singly- 
thermostated oscillator equations are not “ergodic”, so that Gibbs’ statistical mechanics is 
unable to describe the oscillator’s properties. For the next 15 years, which included many 
failed attempts, no singly-thermostated oscillator models were found to be ergodic. 

This letter announces our recent achievements toward the longstanding goal of ergodic 
singly-thermostated oscillator models. We have carried out a comprehensive exploration of 
a previous model claimed to be ergodic, and found that it is not. As a result of those inves¬ 
tigations we have found a path leading to a singly-thermostated and physically motivated 
ergodic model for the harmonic oscillator. We lay out the details of these discoveries in what 
follows and encourage the reader to help explore the new areas opened up by our work. 


II. ERGODICITY IS TYPICALLY ABSENT IN THE SF MODEL 

In 2001 Sergi and Ferrario [ SF ] announced that they had found an ergodic thermostated 
oscillator model^. In addition to the oscillator coordinate, momentum, and thermostat 
variable {q,p, () their model includes a parameter ly which can be either positive or negative: 

Q = P(^ + ; P = -Q - (P C = [P^ -T -qpp ]/t^ ] V = C ■ 

Here, and in what follows, we will ignore the fact that SF actually solve the above four 
equations, not just the three shown below: 

{ q = p{l + (u) ; p= -q-Cp; t = P^ - T - qpu } [ SF ] , 
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This is because their work was based on a Hamiltonian with two degrees of freedom. Consider 
a particular initial condition {q,p,C) that evolves in some time f to a unique {q',p',C)- 
The latter variables do not depend on the initial value of p, which could be given or not, 
arbitrarily. The fourth equation, for the evolution of a variable which is the time integral 
of ( , plays no role at all in the dynamics of {qpQ and can so be ignored, which we do 
throughout. This extraneous variable obscured the fact that SF implicitly claimed ergodicity 
for a sm^'/^thermostated oscillator. As a result, this desirable feature of their relatively 
widely-cited paper has been previously ignored. However, as a consequence of removing p, 
the symplecticity of the dynamics disappears. 

Like the NH model the SF oscillator has mirror symmetry {+q, +p) <—> (—g, —p) . In 
addition the time reversibility of the Sergi-Ferrario equations requires that the functions 
p and ( , as well as the parameter v , all change sign in the reversed motion with the 
coordinate values unchanged. For clarity we have replaced Sergi and Ferrario’s parameter 
“r” by V throughout the present work. This change emphasizes that an increase in | z/ | 
reduces the response time of the thermostat terms. 

For the remainder of this study, we choose to keep r = 1. Usually r , which represents the 
relaxation time of the dynamics, is chosen according to the relatioi>I^: = kTjuP'^ where 

bj is the angular frequency of the system. In our present case, since the system comprises 
a single harmonic oscillator with unit mass and spring constant, a; = 1. Additionally, most 
of the work ascertaining the ergodicity of thermostatted dynamics has taken the relaxation 
time to be unity. We wish to highlight the fact that if the relaxation time is chosen too large, 
it will have no effect on the system dynamics, while if r is chosen too small, the equations 
become too stiff. 

Sergi and Ferrario claimed that their four [ but actually only three, for the reason just 
cited ] oscillator equations^ were ergodic ( hlling out the entire three-dimensional Gaussian 
distribution ) for v > 0.5 . That surprising claim sparked the present work. To begin 
our exploration of their model we carried out a simulation of the SF equations with the 
temperature T and parameter v both equal to unity and with the initial conditions (g, p, C) = 
(1,1,1) . Figure 1 shows the resulting torus, colored according to the local flow instability. 
Evidently this special case of the SF model is definitely not ergodic. 

The difficulty in isolating a small embedded torus by looking at the global dynamical^ 
prompted us to investigate the Poincare section at C = 0. In fact, any other typical Poincare 
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Figure 1: The torus shown here results from the intial conditions {q,p,C) = (Ij 1) 1) using Sergi 
and Ferrario’s original equations with ly = +1 . The local values of the largest Lyapunov exponent 
on the torus are indicated by color : —1.06 (blue) < Ai(t) < 1.89 (red) . Its time-averaged mean 
value, Ai = ( Ai(t) ) is zero. The temperature is unity. 

section would have served our purpose. Recall that Gibbs’ probability density is Gaussian 
in both q and p. Accordingly sections in q and p (as well as in (C) that are far from origin 
are atypical, and may not give any useful results. So long as the section chosen is a typical 
one, the dynamics within it can be studied to understand ergodicity. 

Rather than abandoning the SF approach we also looked for modifications that might 
be ergodic. Ghanging the parameter u from 1 to 2 or 3 or 4 or 5 or 6 and applying due 
diligence led in each case to the discovery of nested tori. Typically the tori penetrate the 
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plane C = 0 in fonr widely-separated distinct places. Figure 2 illustrates these “period- 
four” equilibrium points for the SF equations. Just as in the other Figures the online version 
is colored according to the local value of the largest of the three Lyapunov exponents, Xi(t). 
We denote the long-time average value of this exponent by Ai = ( Ai(t) ) . 

Holes in the chaotic sea are most easily found visually. Then, zooming in on such a hole 
the central point corresponding to a periodic orbit can be found. By first looking at cross 
sections decorated by a million penetration points and then zooming in on the holes we can 
obtain precise six-hgure estimates for the (g, p, 0) points that lie at the center of each hole, 
on the central periodic orbit. Viewed in the {q,p, 0) plane, diligent searches showed that the 
six choices of u shown in Figure 2 , all include simple tori centered on a periodic orbit and 
embedded in a chaotic sea. Looking at the Figure, the relatively small but clearly visible 
holes can be seen for z/ = 2,4, 5, and 6 . The large irregular holes for u = 1 form a cross 
section of the torus shown in Figure 1. The four tiny holes corresponding to = 3 are too 
small to see without zooming in. Some of the details of these investigations are described in 
what follows, along with a concluding Summary, Discovery, and Advice section. 

III. LYAPUNOV INSTABILITY AND GAUSSIAN MOMENTS 

The largest time-averaged Lyapunov exponent Ai measures the exponential tendency for 
two nearby chaotic trajectories to separate, S{t) ~ (5(0)e+'^* . The local value Ai(t) exhibits 
fluctuations, even in the regular toroidal regions where the long-time average, Ai , is zero. Ai 
is positive in the chaotic sea. It is a measure of the chaos there. For online viewing we have 
included the local values with color, ranging from blue to red as the exponent increases. For 
the model of Figure 1 the long-time-averaged exponent is equal to zero, as expected for a 
two-dimensional torus in a three-dimensional space. Simulations with z/ = 2, 3,4, 5, 6 looked 
much more promising, as they all generated “fuzzy balls” in (g,p, C) space. We picture the 
three-dimensional Gaussian distribution, proportional to ^ as a fuzzy ball. It 

is evident that the density falls off exponentially in all directions as one moves away from 
the “center” of the pictured ball. 

We next investigated the ergodicity of these fuzzy balls by measuring the time-averaged 
moments { ( g^, g^, ) } . For every value of v , using a spacing of 0.05 with 

0 < z/ < 6 , we found that the deviations from the even Gaussian moments are small 
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Figure 2: Cross sections at C = 0, for different values of u, with the local values of the largest 
Lyapunov exponent Ai(t) colored from ~ —1.0 (blue) to ~ +1.0 (red). The top, middle and 
bottom figures on the left correspond to i/ = 1,3 and 5, respectively. Likewise, the figures on the 
right correspond to nu = 2,4 and 6, respectively. The white curves correspond to intersections 
with the nullcline surface. There the phase-point velocity is parallel to the {q,p,0) plane with 
q = — 1)/ i^p) ■ The tori for z/ = 2,4, 5 and 6 penetrate the plane within four roughly-triangular 

holes in those cross sections. The tori for u = +3 are much smaller, as is detailed in Figure 3. The 
four penetrations occur as two pairs of mirror-image points : {q,p,0) = (±1.04031, ±0.39432, 0) 
and (±1.18488,+0.94527,0) . The temperature is unity throughout. 
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and masked by fluctuations whenever the tori diameters are small. These deviations led us 
to a topological investigation of the distributions f{q,p,C) different cross sections, that 
constitutes a much more sensitive indicator of nonergodicity than either the moments or 
the one-dmensional probabilities associated with them. We carried out visual inspections 
of zero-C cross sections like those shown in Figure 2. For values of u both smaller and 
larger than the borderline value u = 0.5 put forward by Sergi and Ferrario, the sections 
all revealed well-defined “holes”. The “holes” that can be seen in the Figure correspond to 
toroidal solutions which penetrate the surrounding chaotic sea. 

Figures 3 and 4, corresponding to u = -|-3 and u = +2.9 respectively, reveal triangular 
regions enclosing nested tori. These tori are a clear proof of nonergodicity. The reversibility 
of the motion equations suggests that changing the signs of iy^pX) the initial condi¬ 
tions will simply reverse the trajectories. Numerical work, using fourth-order, fifth-order, 
and adaptive Runge-Kutta integrators bears that expectation out. The rotation of the tri¬ 
angle seen in the closeups ( with linear zooms of factors of ten and one hundred ) from 
longest-side-“up” to longest-side-“down” suggests a singular region in between, which fur¬ 
ther investigation locates near v = -1-2.903521 . 

It appears that the tori shrink to a single periodic orbit at this value of v before enlarging 
again as v increases further. The periodic orbit is shown in Figure 5. A zoom into this 
region by a factor of ten million places an upper limit of 3 x 10“® on the size of the thin torus 
that presumably surrounds the periodic orbit. The limiting torus has a winding number of 
(1/3), which means that the orbit rotates through an angle (27r/3) the short way around 
the torus for each time around the long way. Since the periodic orbit is a neutrally stable 
fixed point in the Poincare section, it is surrounded by a region that very slowly Alls in by 
orbits approaching from the chaotic sea that spend a long time in its vicinity, an example 
of which is evident in Figure 4. Thus it appears that at this singular value of v the system 
may be ergodic but only after an infinite time. 

It seems likely that the mathematical [ as opposed to physical ] form of the oscillator 
equations, where p differs from q and where the friction coefficient depends on qp , was 
offputting for later investigators so that these 2001 contradictions with the literature of the 
1990s passed either unnoticed or at least undeclared until now. 

The systematic explorations carried out by Bauer, Bulgac, Ju, and Kusnezov, for simple 
systems including the harmonic oscillator, suggested that quartic thermostating terms like 



Figure 3: An enlarged version of the v = +3 case shown in Figure 2 shows four barely-visible 
“holes” in the C = 0 cross section. Magnification of the lower-left hole, tenfold in both the q and 
the p directions, reveals a sharp triangular boundary between nested tori on the inside and chaos 
on the outside. Because the location of the toroidal region is a smooth function of v successive 
approximations track its center to z/ = -|-2.9 , shown in Figure 4, and finally to u = -1-2.903521 , 
where the sidelength of the triangle is less than 3 x 10“® . Color indicates the local value of the 
Lyapunov exponent Ai(t) , with red positive and blue negative. 

—(p^ or —(^p best promote chaos^^— . Accordingly, we modified the Sergi and Ferrario 
eqnations to inclnde a cnbic ( rather than linear ) dependence on the friction coefficient : 

{q = p{l + ; p= -q-C^p] ( = p'^ - T - qpu } [ SF' ] . 
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Figure 4; A hundred-fold linear zoom in q and p reveals the toroidal solutions within a triangular 
reqion of width of order 0.001 . The white lines indicate {q,p) tracks of points moving parallel to 
the plane. Red and blue represent positive and negative local Lyapunov exponents Xi{t) . The 
cross section displays the mirror ( or inversion ) symmetry {+q, +p) <—>■ {—Q, —p) ■ is +2.90 
here. 

We also followed References 8-10 by considering a quadratic version with —)■ ^|(^| . The 
reader can check to confirm that these equations for the flow velocity still satisfy the sta¬ 
tionary phase-space flow equation , (df/dt) = 0 = —V • {fv) . Here v = {q,p,C) is the 
three-dimensional flow velocity in (qpC) space. The probability densities for the quadratic 
and cubic generalized friction coefficients C vary as ^ and . 
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Figure 5: The stable periodic orbit near v = 2.903521 and its projections into the {q,p), iq,C)^ 
and (p, C) planes. This orbit has a period of 12.2945528 and crosses the C = 0 plane near {q,p) = 
(+1.04099057, -0.395775939) . 

IV. TWO-THERMOSTAT HARMONIC OSCILLATOR MODELS 

Beginning about 1990 two-thermostat models were developed, mostly based on controlling 
pairs of moments^^— . Applications established the mathematical consistency of such models 
with Gibbs’ canonical distributions, with barrier-crossing problems, and with Brownian 
motion problems. The Ju-Bulgac, Martyna-Klein-Tuckerman, and Hoover-Holian models 
all generated the canonical distribution. The controversial ergodicity of the MKT oscillator 
was investigated, and conhrmed with particular care, as the 2014 Snook Prize problem^^i^. 
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Rather than formulating two controls over the potential or kinetic energy this MKT “chain” 
model used a second thermostat variable ^ to thermalize the fluctuations of the first, : 

{q=p; P = -q-Cp; C = p"-T-^C; ^ = C" - T } [ MKT ] 

In 2014 Patra and Bhattacharya discovered that the doubly-thermostated oscillator equa¬ 
tions , 

{q = p-^q; p=-q-Cp; C = p"-T; i = q^-T}[FB = SE], 

were not ergodio^^. By coincidence Sergi and Ezra had already found this same result in 
20101^. We discovered this by noticing that their Figure 2 looked identical to Patra and 
Bhattacharya’s Figures 2c and 2d in Reference 15. The key to understanding ergodicity and 
its lack in these simple oscillator systems lies in distinguishing two qualitatively different 
types of “holes” in the cross sections of the flow. We turn to that next. 


V. HOLES IN THE SINGLY-THERMOSTATED CROSS SECTIONS 

The holes found here in the cross sections are reminiscent of those found recently by Patra 
and Bhattacharya!^. They investigated the two unstable fixed points generated by the four¬ 
dimensional MKT oscillator equations. Evidently the centers of the largest holes found in 
the present work typically include fixed cycles of four repeating points of the mapping from 
one penetration of the plane at C = 0 to the next (three intermediate penetrations separate 
pairs of point repetitions ) . The holes are especially clear for the case z/ = 2 in Figure 2. 

Figure 4 illustrates a relatively sensitive case, u = -1-2.90 , using the original Sergi- 
Ferrario equations with linear friction. Apart from four tiny similar holes in the section, the 
chaotic sea outside them has a Gaussian distribution. This is consistent with the largest of 
the long-time-averaged Lyapunov exponents ( as well as with the complete spectrum of four 
exponents ) vanishing for all those trajectories which pass through the holes. 

We took the precaution of solving this problem with three different integrators ( fourth- 
order, hfth-order, and adaptive Runge-Kutta ) and a variety of fixed and variable timesteps, 
all in a diligent effort to avoid numerical errors. For a purely-Hamiltonian harmonic oscillator 
it is well-known that the fourth-order method gradually loses energy while the fifth-order 
method gains. The good agreement of all three integrators with one another shows that the 
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Figure 6: Symmetry breaking at = ±2.9 . The time-reversible (g, ±p, ±C, ±z^) trajectories 
trace out identical coordinate sequences, but in the opposite time direction. Even so, differences 
in the recent past histories of the forward and reversed trajectories lead to the totally different 
local values of the largest Lyapunov exponent Ai(t) shown here. The global forward and backward 
time averages of the local exponents Ai match. The section for u = —2.9 has been reflected about 
the coordinate axis {p = 0) to show that the trajectory penetrations as well as the white nullcline 
intersections parallel to the ((7,p, 0) plane are identical in the forward and reversed trajectories. 

nonlinearities of the differential equations dominate the errors ( on the order of 10“^^ or less 
at each timestep ) from the finite precision of the simulations. By simply searching for holes 
and evaluating the largest Lyapunov exponents within them, or by evaluating the largest 
Lyapunov exponents for millions of randomly-chosen initial conditions it is relatively easy 
to separate the chaotic and quasiperiodic regions. 
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VI. SUMMARY, DISCOVERY, AND ADVICE 


Three families of singly-thermostated oscillators, with linear, quadratic, and cubic fric¬ 
tion, were formulated to obey Gibbs’ canonical distribution for an oscillator. All provided 
chaos but none was ergodic. In the linear case it is possible that there is a ergodic solution 
in the vicinity of = 2.903521 . We can place an upper limit, ~ 10“^® , on the nonergodic 
measure there. In every case we examined ( with independent calculations in India, Nevada, 
and Wisconsin ) the cross sections revealed mixed solutions, holes containing nested tori 
embedded in a chaotic Gaussian sea. Figure 5 shows such a stable torus. 

This situation bears a qualitative resemblance to solutions of the original NosGHoover 
oscillator. Visualization and the local values of the largest Lyapunov exponent are the two 
most valuable tools for distinguishing the two solution types. To evaluate the local exponent 
Ai(t) requires both a “reference” trajectory and a nearby “satellite”. This increases the 
computer time required by only a factor of two. It is convenient to rescale the separation 
between the two similar trajectories ( by displacing the satellite toward the reference 
at every timestep so as to determine Ai(t) . 

These three-dimensional oscillator problems help illuminate features of ergodicity searches 
for the four-dimensional flows obtained with two thermostat variables. Our work here has 
shown that the Sergi-Ferrario oscillators are at best seldom ergodic ( assuming only that the 
hundreds of cases we examined are typical ). 

A particularly fascinating aspect of the fully time-reversible Sergi-Ferrario model is the 
symmetry breaking illustrated in Figure 6. The forward structure of the flow’s Lyapunov 
instability is much more complex than the totally different backward structure despite the 
fact that any trajectory obeying the SF motion equations can be followed just as well 
backward as forward. This “Arrow of Time” asymmetry of the largest Lyapunov exponent 
Ai(t) for a simple fully time-reversible flow deserves further study. 

We recommend the challenge of taking up the search for ergodic three-dimensional oscil¬ 
lator models. In pursuing this elusive goal it seems to us highly desirable to maintain the 
conventional relation between the coordinate and the velocity, q = p ■ It is also desirable to 
resist such physically-artihcial accelerations as the qp contribution to the thermostat vari¬ 
able ( . At their best control variables should utilize transparent and meaningful origins. It 
is certainly possible that with increasing computer power searches such as those carried out 
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by Sprotlji^ could uncover a host of new variations of the Sergi-Ferrario or Sergi-Ezra-Patra- 
Bhattacharya equations which are simultaneously robust, useful, physically meaningful, and, 
above all, ergodic. 

In closing, we have recently discovered a particularly promising direction embodying 
“weak control” of the momentum through a choice of the parameters (a,/ 3 , 7 ) in the set of 
three moment-based equations of motion : 

{i = +p \ p = -9 - Cl ap + /^(p“/r) + 7(p^/r^) 1; 

C = a| (pVt) -1 ] + /3[ (pVt") - 3(p7r) ] + 7[ (p«/r7 - 5(pVr7 ]}. 
Computational searches in (a, /5, 7 ) space suggest that there are regions where the singly- 
thermostated oscillator samples the entire Gibbs’ distribution. One such combination which 
appears to be ergodic is {a,{3,'y) = (1.50,0.00,-0.50) . We suspect there are many more. 
Finding them will conclude a search set in motion by Shuichi Nose some thirty years ago. 
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